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ALGORITHMS  FOR  SINGLE-SIGNAL  AND  MULTISIGNAL 
MINIMUM-CROSS-ENTROPY  SPECTRUM  ANALYSIS 


INTRODUCTION 

Multisignal  minimum-cross-entropy  spectral  analysis  (multisignal  MCESA)  is  a  method  for 
estimating  the  power  spectrum  of  one  or  more  independent  signals  when  a  prior  estimate  for  each  is 
available  and  new  information  is  obtained  in  the  form  of  values  of  the  autocorrelation  function  of  their 
sum  (1  ].  This  method  is  a  generalization  of  MCESA  [2]  and  reduces  to  it  when  there  is  only  one  sig¬ 
nal.  One  application  of  multisignal  MCESA  is  io  noise  reduction  [3],  If  one  obtains  autocorrelation 
measurements  for  a  signal  with  independent  additive  interference,  and  if  one  has  some  prior 
knowledge,  expressible  as  spectral  estimates,  concerning  the  signal  and  the  noise,  the  method  yields 
new  signal-  and  noise-spectrum  estimates  that  take  both  the  prior  estimates  and  the  autocorrelation 
information  into  account. 

We  here  present  two  algorithms  that  implement  multisignal  MCESA.  One  is  slow,  but  general; 
the  other  is  considerably  faster  and  applies  to  an  important  special  case. 

The  first  algorithm  accepts  as  inputs  prior  estimates  of  the  power  spectra  of  one  or  several 
independent  signals  and  measured  autocorrelation  values  for  the  sum  of  the  signals,  and  it  produces  as 
outputs  posterior  spectrum  estimates  that  are  consistent  with  the  autocorrelation  estimates.  The  spectra 
are  represented  by  discrete-frequency  approximations  —  lists  of  spectral  powers  at  specified  frequencies. 
The  algorithm  implements  the  methods  of  Refs.  1  and  2  in  full  generality;  the  discrete  frequencies  may 
be  arbitrarily  spaced,  and  the  autocorrelations  may  be  given  at  arbitrarily  spaced  lags.  The  algorithm  is 
iterative;  at  each  iteration,  it  computes  frequency-domain  quantities  (spectra)  and  time-domain  quanti¬ 
ties  (autocorrelations  and  Lagrange  multipliers)  that  are  related  to  each  other  by  procedures  equivalent 
in  complexity  to  the  naive  (slow)  algorithm  for  the  discrete  Fourier  transform. 

The  second  algorithm  applies  to  the  special  case  in  which  the  prior  spectral  estimates  are  of  the 
all-pole  form  that  results  from  maximum-entropy  spectra!  analysis  (MESA)  [4,5]  or  linear-predictive 
coding  (LPC)  [6].  Such  spectra  may  be  represented  in  various  ways  by  finite  families  of  parameters 
that  determine  them  uniquely  —  for  example  by  a  gain  and  inverse  filter  coefficients  or  by  a  gain  and 
reflection  coefficients.  In  particular,  the  spectra  may  be  represented  by  a  finite  number  of  autocorrela¬ 
tion  values.  This  algorithm  accepts  prior  autocorrelation  estimates  for  the  individual  signals  in  place  of 
prior  spectral  estimates.  The  prior  autocorrelation  estimates  for  the  individual  signals  and  the  measured 
autocorrelations  of  the  sum  of  the  signals  are  given  at  equispaced  lags  beginning  at  zero.  The  posterior 
spectral  estimates  that  result  from  data  of  this  form  are  all-pole  spectra  like  the  prior  estimates,  though 
not  necessarily  of  the  same  order.  The  outputs  from  the  algorithm  may  be  posterior  autocorrelation 
estimates  or,  alternatively,  any  of  various  equivalent  families  of  LPC  parameters. 

The  next  section  gives  enough  background  for  a  statement  of  the  problem  to  be  solved;  for  a 
fuller  discussion,  see  Refs.  1  and  2.  The  third  section  presents  the  first  algorithm;  the  fourth  section 
presents  the  second  algorithm.  The  fifth  section  contains  a  general  discussion. 


Manuscript  approved  November  2.  1982. 
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THE  PROBLEM 

Consider  a  number  p  of  unknown  real,  independent,  band-limited,  discrete-spectrum  signals  with 

no  dc  component;  let  their  spectra  be  Sh>  h  =  1 . p.  with  spectral  powers  Shk  at  frequencies  ±,f\. 

where  /*  >  0  for  it  =  1,  .  .  .  ,  n.  Let  Phk  be  a  prior  estimate  for  Shk.  Now  suppose  we  are  given  values 
R,[0{  of  the  autocorrelation  function  of  the  sum  of  the  signals  for  lags  r,,  r  =  0,  .  .  .  ,  m.  Under  the 
assumption  that  the  processes  are  independent,  we  write  Rru"  as  a  sum  of  individual  autocorrelations, 

Rr'°'=  ±R„, 

h~\ 

where 


Rfir  =  £2SM.  cos  2nfktr. 

k- 1 

We  wish  to  use  the  autocorrelation  information  to  obtain  improved,  posterior  estimates  Qh  of  the  Sh. 
The  multisignal  minimum-cross-entropy  estimate,  given  by  Eq.  (18)  of  Ref.  1,  is 


Qhk  -  Qkkifi)  =  QkM i.  -  •  •  .fiJ 

l 


l  m 

-jr-  +  I  PrCrk 

rhk  r-0 


(1) 


where 


Crk  =  2  cos  2nfktr.  (2) 

The  parameters  /3r,  r  —  1,  .  .  .  ,  m,  are  Lagrange  multipliers  associated  with  a  constrained  minimization 
problem;  the  cross  entropy  of  a  pair  of  probability  distributions  related  to  the  Ph  and  the  Qh  is  minim¬ 
ized  subject  to  constraints 

Rr'°l=  t  t&kWCrk  (3) 

h- 1  *-l 

( cf.  Ref.  1,  Eq.  (21)]  that  require  the  posterior  estimates  Qh  to  be  consistent  with  the  given  values  of 
the  summed  autocorrelations.  The  /3r  are  to  be  chosen  so  that  Eq.  (3)  is  satisfied.  The  first  algorithm 
applies  in  this  case,  accepting  the  quantities  tr ,  Rrl0\  f, ,  and  Phk  as  inputs  and  computing  the  Qhk  as 
outputs. 

In  the  case  to  which  the  second  algorithm  applies,  we  take  the  spacing  between  autocorrelation 
lags  as  a  unit  of  time  and  its  reciprocal  as  a  unit  of  frequency.  We  thus  write  tr  =  r,  r  =  0,  .  .  .  ,  m.  In 
place  of  a  discrete-frequency  approximation,  we  write  spectra  as  functions  of  a  continuous  variable  /  in 
the  interval  <  /  <  'h.  Thus  we  have  prior  estimates  Pk(f)  of  the  spectral  power  densities  Sh(f). 
Corresponding  to  Eq.  (1)  we  have 

<?.</> - j - i - ■  <4' 

p  ,  77  +  2  £  Pr  cos  2 nfr 

rhyJ 9  r-0 

and  the  constraints  to  be  satisfied  are 

Rr"n  =  2  <?„(/)  cos  2n.fr  df. 

The  second  algorithm  avoids  direct  computation  with  spectra  Ph  and  Qh  —  that  is,  it  avoids  the 
use  of  discrete-frequency  approximations  to  Ph  and  Qh  and  numerical  integration  over  frequency.  This 
is  possible  since,  as  we  will  see,  all  the  spectra  that  occur  are  of  the  all-pole  (MESA  or  LPC)  form  [5,6] 
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<?*(/)  = 


Eh 

T,  ~  p" 

e2"1" 

1-0 


where  the  a/,,  are  inverse  filter  coefficients  and  Eh  is  a  gain.  This  is  equivalent  to 

(?*(/)  =  — - 1 - • 

2^  kf,,  cos  2rr /r 
(-0 


(5) 


(6) 


where  the  Lagrange  multipliers  \h,  are  related  to  Eh  and  the  ahl  by 

1  M 

^  HO  r-  ^  .  @hs 

2Eh  sto 

and 

j  M-t 

^  hi  r  ^E^hs^h.s+t1  t  1 . 

Lh  5—0 


(7) 


The  prior  spectra  Ph  have  the  all-pole  form  by  assumption;  we  express  them  in  the  form 


n</>  -  — — 1 - 

2X  ** cos  lnf! 


We  may  assume  that  the  prior  spectra  are  of  order  M  ^  m ,  for  we  can  extend  the  sequences 
\h0,  XM,  .  .  .  with  zeros  if  necessary.  To  see  that  the  posterior  spectra  must  have  the  same  form,  we 
substitute  this  expression  for  Ph  in  Eq.  (4).  We  find  that  Qh  is  given  by  Eq.  (6)  with 


[^)w* 


0  ^  t  <  m, 
m  <  t  ^  M. 


(8) 


The  autocorrelations  Rh,  of  Qh  are  related  to  the  \  h,  by 

Rh,  =  f0h2Qh{f)  cos  2n ft  df  (9) 

f  'h  cos  2-rrft 

Jo  M  aJ' 

£ \hu  cos  2t t/u 

u—0 

and  the  constraints  can  be  expressed  as 

Y,Rhr=Erw,  r  —  0 . m.  (10) 

h 

Either  family  of  parameters,  Rh  =  (Rh o,  Rh\ . R or  kh  =  (kh0,  kh\ . khM),  is  sufficient  to 

determine  Qh ,  and  there  are  algorithms  to  compute  either  family  from  the  other  without  numerical 
integration  of  the  right  side  of  Eq.  (9).  Thus  we  can  work  in  terms  of  Rh  and  A*  rather  than  Qh. 

The  second  algorithm  begins  with  the  summed  autocorrelations  R,'nl  and  the  autocorrelations 

Rh,  =  2 JQ  Ph(f)  cos  2nft  dt 

of  the  prior  spectral  estimates  Ph.  It  computes  the  autocorrelations  Rh,  of  the  posterior  spectral  esti¬ 
mates  Qh ;  the  Lagrange  multipliers  \h,  of  Eqs.  (6)  and  (9)  are  also  available  as  outputs. 
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THE  FIRST  ALGORITHM 
Derivation  of  the  First  Algorithm 

When  Eq.  (1)  holds,  the  right-hand  side  of  Eq.  (3)  for  each  r  is  a  function  of  the  Lagrange  multi¬ 
pliers. 

w-ti-, — ^ — c'<-  <n» 

4-  +  z 

rhk  j-0 

we  need  to  solve 


Frip)  =  Rr10',  r  =  1 . m. 


(12) 


for  p.  The  Newton-Raphson  method  for  approximating  solutions  to  such  systems  of  equations  starts 
with  an  initial  guess  p{0)  at  the  solution,  iteratively  computes  approximations  /J(I),  0(2) . by 

0<'i  =  p<'-D  -i-  &p(i),  (13) 

where  A  PU)  satisfies 

m  BF(B{i~'l>) 

F,(pu~")  +  T  -I,—-  A p™  =  Rrl0\  (14) 

st f  dpl'~u 

and  stops  when  the  p(l>  have  satisfied  some  convergence  criterion.  In  matrix  notation,  for  fixed  i  the 
solution  W{>)  of  Eq.  (14)  is  given  by 

A  P{i)  =  M-1  v, 


where  v  is  the  vector  with  components 

and  M  is  the  matrix  with  elements 


v,  =  Rrw  -  Fr(p(i-l)), 
BFr(pu~") 


= 


9/3 


(<-D  ' 


Differentiating  Eq.  (11),  we  obtain 


p  n 

I  I 

A— I  *-t 


Crk  Csk 


m 


—  +  '£tp;l~"clk 


hk  i-0 
P 


I  -crkcsk  x  Qkk(fiu-"y 

*-l  h- I 


(15) 

(16) 


It  follows  that  M  is  given  by 


M  =  -NN\ 


where  the  prime  denotes  transposition  and 

Nrk  =  Crk 


I  Q*<F‘-U)‘ 

h- 1 


(17) 


(18) 
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From  Eqs.  (11),  (1),  and  (18)  we  obtain 


F'W-")  -  I-V.a  , 

A-'  i  e**(0,,'-1,‘2 


We  may  therefore  express  the  vector  v  defined  in  Eq.  (16)  by 

v  =  -  N  w, 

where 


~-v<  + 


I(?W(/9,"n)2 


and  the  x*  are  any  solution  of 


(19) 


(20) 


(21) 


xk 


10«0<'-u)2 


=  /?;- 


By  Eq.  (18),  this  is  equivalent  to 


£  Qx,  =  /?;■ 


or,  in  matrix  notation, 


Cx  =  tf10'. 


(771 


Substituting  Eqs.  (17)  and  (20)  in  Eq.  (15),  we  obtain 

A0(,)  =  (NN')-'Nw. 

The  quantity  (NN’)-1N  (when  NN'  is  nonsingular)  equals  the  Moore-Penrose  generalized  inverse  [7] 
N’^ofN'.  Thus, 

A  fi,n  =  N'V  (23) 

(This  algorithm  was  first  coded  in  the  programming  language  APL.  Equation  (23)  is  in  a  convenient 
form  for  translation  into  APL,  since  there  is  an  APL  primitive  function  [81  that  will  yield  the  general¬ 
ized  inverse  of  a  matrix.)  We  can  express  the  solution  of  Eq  (22)  in  terms  of  the  generalized  inverse. 
We  assume  that  Eq.  (3)  provides  m  independent  conditions  on  the  Qh.  thus  the  rows  of  C  are  linearly 
independent.  This  implies  that  the  generalized  inverse  is  a  right  inverse;  then  CC'  is  an  identity  ma¬ 
trix.  It  therefore  suffices  to  set 

x  =  C7T0'  (24) 


for  a  solution  of  Eq.  (22)  The  essential  ingredients  of  the  algorithm  are  now  at  hand.  We  choose  an 
arbitrary  starting  value  for  f}10’,  compute  the  by  Eq.  (1).  and  define  x  by  Eq.  (24>.  We  then 

enter  a  loop  in  which,  on  the  /th  iteration,  we  first  compute  p(l>  by  Eqs.  (13),  (18),  <21 ),  and  (23)  and 
then  compute  by  Eq.  (1).  The  loop  terminates  when  the  have  converged  by  some 

appropriate  criterion.  In  practice,  some  modifications  of  the  procedure  are  necessary  to  prevent  the 
Newton-Raphson  algorithm  from  behaving  badly  in  case  the  arbitrary  starting  point  is  not  near  enough 
to  a  solution.  These  are  incorporated  in  the  following  step-by-step  outline  of  the  algorithm  as  the 
adjustment  to  A  P  in  Step  6. 
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Summary  of  First  Algorithm 

In  what  follows  we  will  use  a  left  arrow  to  denote  assignment  to  a  variable  of  a  new  value,  possi¬ 
bly  computed  in  terms  of  the  current  value.  Thus,  corresponding  to  Eq.  (13)  we  write  p<—p  +  \p 
(without  superscripts)  in  Step  7.  The  inputs  to  the  first  algorithm  are  three  vectors,  r,  R'°',  and  /,  and 
a  matrix  P.  listing  respectively  the  lags,  autocorrelations,  frequencies,  and  prior  spectral  estimates,  tr, 
Rr'°\  fk,  and  Phk,  r=  1 . m.  k=  1 . n\  h=\,  ...  ,p.  The  output  is  a  matrix  Q  listing  the  poste¬ 

rior  spectral  estimates  Qhk. 

The  steps  of  the  algorithm  are  as  follows: 

Step  1.  Compute  the  matrix  C  by  Eq.  (2)  and  compute  the  vector  x  by  Eq.  (24). 

Step  2.  Assign  initial  values  to  variables 

p  -  (0 . 0) 

Q-P. 

Step  3.  (Begin  main  loop.)  Save  the  current  value  of  Q  for  later  comparison  with  the  new 
value: 

Q°  Q- 

Step  4.  Compute  a  tentative  value  of  A/J  by  Eqs.  (18),  (21),  and  (23). 

Step  5.  Compute  tentative  values  Thk  for  1  /  Qhk  by 

Thk  —  +  A/?,)cy* 

r hk  r 

[cf  Eq.  (1)]. 

Step  6.  If  Thk  >  0  for  all  h  and  k ,  proceed  to  Step  7.  Otherwise,  adjust  A/B  by  the  replace¬ 
ments 

0.9A/3, 

A/T  -  - .  n0T 

1  -  min  QSc  Thk 

hk 

and  recompute  the  Thk  by  Step  5  before  going  on  to  Step  7.  (The  denominator  on  the 
right  side  was  chosen  to  make  the  recomputed  values  of  Thk  nonnegative.  The  factor 
0.9  is  included  to  assure  strict  positivity.  The  value  0.9  is  rather  arbitrary;  in  principle 
it  could  be  any  positive  number  less  than  1.) 

Step  7.  Assign  new  values  to  variables 

P'—p  +  bfi 

and 

Qhk  1/  Thk  ■ 

Step  8.  If  the  new  value  of  Q  equals  the  previous  value  (saved  as  Q°)  to  within  a  specified 
relative  tolerance,  terminate  the  algorithm  and  return  Q  as  the  result.  Otherwise  go 
back  to  Step  3  and  repeat  from  there. 
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THE  SECOND  ALGORITHM 

In  this  subsection  we  sketch  a  multisignal  MCESA  algorithm  that  works  in  terms  of  autocorrela¬ 
tions  and  Lagrange  multipliers.  The  next  three  subsections  discuss  useful  standard  subroutines  and 
derive  in  detail  the  major  steps  of  the  algorithm.  The  final  subsection  presents  the  algorithm  in  its 
entirety. 


The  algorithm  begins  with  the  measured  autocorrelations  R'J'  of  the  summed  signals  and  auto¬ 
correlations  /?*,  of  the  prior  spectral  estimates  Ph.  It  is  assumed  that  the  R'01  are  given  at  equispaced 

lags  r  =  0.  1 . nr,  the  Ph  are  all-pole  spectra;  and  the  Rj,\  are  given  at  equispaced  lags 

t  =  0,  I . Mh,  where  Mh  is  at  least  as  great  as  the  order  of  the  all-pole  spectrum  Ph.  There  is  no 

loss  of  generality  in  assuming  that  Rfi  is  given  out  to  lag  M  >  m  for  all  h ,  since  then  are  standard 
LPC  methods  for  extrapolating  autocorrelation  sequences  —  see  the  following  subsectk  .  1  the  sub¬ 
section  "Summary  of  Second  Algorithm."  From  the  families  R®  of  prior  autocorrelatio  ?£,  one  can 
obtain  the  corresponding  Lagrange  multipliers  A®,;  these  satisfy 


cos  2n.fu 


u- 0 


cos  2rr  J'u 


df 


[cf  Eq.  (9)].  Initial  trial  values  are  given  to  the  posterior  autocorrelations  Rhl,  the  corresponding 
Lagrange  multipliers  khl,  and  the  parameters  @r.  {Initial  value  assignments  Rh  —  R®,  kh*~ 

R  —  (0 . 0)  are  reasonable  and  appear  to  work  well  in  practice.)  These  values  are  then  adjusted 

iteratively  with  the  aim  of  simultaneously  establishing  Eqs.  (8),  (9),  and  (10).  At  each  iteration.,  the 
Aa  are  recomputed  from  the  current  values  of  Rh,  Eq.  (9),  which  determines  the  functional  depen¬ 
dence  between  Rh  and  \h,  is  thus  maintained.  Adjustments  ±Rh,  aimed  at  establishing  Eq.  (8),  and 
aimed  at  establishing  Eq.  (10),  are  computed  by  a  linear  approximation  involving  the  partial  deriva¬ 
tives  BRh,/dkhu.  The  replacements  Rh  —  Rh  +  ±Rh,  R—RE\R  are  made,  the  kh  are  recomputed, 
and  the  process  is  repeated  with  the  new  trial  values  of  R,r  Xfl.  and  R  The  repetition  continues  until 
the  values  of  Rh  are  in  sufficiently  good  agreement  on  two  successive  iterations. 


Some  LPC  Algorithms 

Certain  standard  I. PC  signal-processing  procedures  can  be  used  as  subroutines  at  several  points  in 
the  algorithm  FORTRAN  implementations  of  all  of  them  arc  incorporated  in  the  subroutine  LPTRN 
in  Section  4.3  of  Ref  9. 


To  obtain  Lagrange  multipliers  from  autocorrelations,  we  can  start  with  the  Levinson  recursion 
[I0i.  This  yields  inverse  filter  coefficients  ah{\  =  I,  ahl .  a,M  and  a  gain  Eh  such  that 


2  Eh  cos  2-nft 

"7T  f dJ 


M 


(25) 


[cf.  Eqs.  (9)  and  (5)].  Then  the  discrete  convolution  of  Eq.  (?)  yields  the  Lagrange  multipliers  khr 
t  —  0 . M.  (The  Levinson  recursion  also  yields  reflection  coefficients  kh\  ,  ...  ,  kMt.  The  condi¬ 
tion  that  | kh! !  <  1  for  all  t  =  1 . M  is  a  useful  criterion  for  checking  that  Rh  is  a  valid  family  of 

autocorrelations— i.e.,  that  it  admits  an  extension  of  positive  type.) 


There  are  inverted  versions  of  the  Levinson  recursion  that  allow  Rh  to  be  computed  from  the 

family  ah  =  (ahn.  ah] . ahU)  of  inverse  filter  coefficients,  either  directly  or  via  the  reflection 

coefficients.  These  will  also  be  of  use. 


Finally,  to  compute  autocorrelation  from  the  Lagrange  multipliers,  we  could  use  the  algorithms 
above  together  with  an  algorithm  for  computing  inverse  filter  coefficients  from  the  Lagrange  multi¬ 
pliers.  There  are  two  procedures  lor  this  —  one  based  on  C  holesky  factorization  of  a  Toeplitz  matrix 
[111  and  one  based  on  Newton's  method  [12.131  We  mention  these  although,  as  it  turns  out.  we  do 
not  actually  need  to  compute  autocorrelations  from  Lagrange  multipliers 
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Partial  Derivative  Calculations 


The  computation  of  A Rh  and  A fi  requires  the  Jacobian  matrices  bRh/dkh  with  elements 
bRhl/bkhu  ( t,u  =  0  ,  . . .  ,  M).  Differentiation  of  Eq.  (9)  yields 


BRh, 

fohu 


J'  —  cos  2  it  ft  cos  2n  fu 

0  M  5 

^khv  cos  2nfv 
i— o 


df 


Since  the  numerator  of  the  integrand  is  equal  to  —  Vi  cos  2n  fit  +  u)  —  xh  cos  2nf\t  -  u|,  we  may 
write 


BRh, 


-s, 


h,t+u 


-  s, 


(26) 


where 


$hv 


cos  2-rrfv 
M 

£X/,,  cos  2ir  ft 
1-0 


df 


v  =  0  ,  ....  2  M.  We  express  this  in  terms  of  the  inverse  filter  coefficients  ahl  and  gain  Eh  correspond¬ 
ing  to  the  Lagrange  multipliers  \hl.  Equation  (7)  [cf.  Eqs.  (5)  and  (6)]  implies  that 


1 


M 


cos  2 nft 

1-0 


2  E„ 
m 

2> i,e2nif' 


It  follows  that 


'h  2El  cos  2 ir/v 


M 
I/-G 


2nift\ 


df 


Regarding  the  ahl  (t  -  0  ,  ....  M)  as  coefficients  of  a  polynomial  A(z)  =  1  +  ah\Z  +  •  •  •  +  aflMzM , 
we  compute  the  coefficients  bhv  of  the  polynomial  A  (z)2  =  1  +  bh]z  +  ■  ■  ■  +  bhi2Mz2M, 


bhv 

l 


(27) 


v  =  0 . 2M.  (The  sum  runs  from  t  =  0  to  v  when  v  ^  M  and  from  t  =  v  —  M  to  M  when 

v  ^  M  )  Then 


v*  2 E2  cos  2 tt  fv 
~2M  F 

I 

w— 0 


df 


Comparing  this  with  Eq.  (25),  we  see  that  the  vector  Sh  depends  on  £*  and  the  vector  bh  in  the  same 
way  that  the  vector  Rh  of  autocorrelations  depends  on  Eh  and  ah.  Thus,  any  algorithm  designed  to 
compute  autocorrelations  from  gain  and  inverse  filter  coefficients  may  also  be  used  to  compute  the  Shv 
from  the  squared  gain  and  the  coefficients  bhv. 


Here,  then,  is  the  procedure  for  computing  bRh/dkh.  Start  with  Eh  and  the  vector  ah  of  inverse 
filter  coefficients;  these  will  be  available  as  a  result  of  their  use  in  the  computation  of  kh  from  Rh. 
Obtain  the  elements  of  bh  by  Eq.  (27).  Compute  Sh  from  Ej2  and  bh  as  just  described.  Finally,  use  Eq. 
(26)  to  obtain  bRh/b\h- 
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Calculation  of  ARh  and  A fi 

We  next  show  how  to  calculate  A Rh  and  A/3  from  the  Jacobian  matrices  dRh/dkh.  We  define  the 
abbreviation  J,  for  the  Jacobians  by 

Jhiu  ~  3  R)„/d\  hu> 

t,u  =  0  ,  ....  M  We  will  also  use  a  rectangular  matrix  D  defined  by 

1 ,  r  =  /. 

Dn  0,  ^  r, 

/•  =  0 . m.  t  =  0 . A/.  In  terms  of  D  and  its  transpose  D',  we  can  express  quantities  such  as 

D Rh,  the  result  of  truncating  the  vector  R  to  its  first  m  +  1  elements;  D'/3,  the  result  of  extending  the 
vector  /3  to  length  M  +  1  by  appending  zeros;  and  DJ,D',  the  (m  +  1)  x  (m  +  1)  submatrix  in  the 
upper  left  corner  of  J,.  We  can  thus  rewrite  Eqs.  (8)  and  (10)  as 

A,  =  A?  +  D'/3 
and 

=  *«>■. 


After  replacement  of  p  and  the  Rh  by  p  +  A/3  and  Rh  +  A Rh,  with  corresponding  replacement  of 
A,,  we  wish  Eqs.  (8)  and  (10)  to  hold,  at  least  to  a  better  approximation  than  before.  Although  kh  and 
Rh  are  nonlinearly  related  by  Eq.  (9),  the  replacement  value  for  A/,,  for  small  changes,  is  given  by 
Aa  +  AAa,  where  AAa  approximately  satisfies 

J.AA,  «  A Rh  (28) 

(That  is,  ^Sio(dRh,/dXhu)^khu  =  A Rhn  t  =  0 . M.)  To  establish  Eqs.  (8)  and  (10)  through  the 

replacement,  we  wish 

kh  +  AA*  =  Aj  +  Dip  +  A/8)  (29) 

and 

D£(/t*  +  A*„)  =  ^,ot  (30) 

h 

to  hold.  By  Eqs.  (28)  and  (29),  A Rh  should  (approximately)  satisfy 

A**  =  Ja(Aa°  +  D'P  -  kh)  +  J*D'A/8.  (31) 

Using  this  equation  to  eliminate  A Rh  from  Eq.  (30),  we  obtain 

D +  Dp, (A?  +  Dp  -  kh)  +  IDJ,D'A/3  =  RM: 

h  h  h 

hence  A/3  may  be  computed  by 

A p  =  |^DJfcD'j-1{/?,m  -  D£[/f„  +  J,(A,°  +  Dp  -  A,)]).  (32) 

The  procedure  for  computing  A p  and  A Rh,  in  summary,  is  first  to  compute  the  quantities 
J,(A°  +  D'/3  -  A,),  then  obtain  A/3  from  Eq.  (32),  and  finally  obtain  the  A/?,  from  Eq.  (31). 

Summary  of  Second  Algorithm 

We  have  discussed  various  parts  of  the  algorithm;  here  is  how  they  fit  together.  Inputs  are  the 
prior  autocorrelation  estimates  for  the  individual  signal  components  h  and  the  measured  autocorrela¬ 
tions  for  their  sum: 


9 


RODNEY  W  JOHNSON 


*/?=<*.'  *A°i . R&4„) 

and 

*ioi  =  (/?r.  *r . /cm>- 

No  initial  assumptions  are  made  about  the  orders  Mh  of  the  prior  estimates. 

The  steps  of  the  second  algorithm  are  as  follows: 

Step  1.  Define 

M  =  max}m,  M\  ,  ....  Mp }. 

Step  2.  By  the  Levinson  recursion,  compute  families  of  inverse  filter  coefficients 
(a/, o=l-  ah j,  .  .  .  ,ahM  )  and  reflection  coefficients  (kh\,  .  .  .  ,khM  )  and  a  gain  Eh 
corresponding  to  R°  for  each  h.  (The  conditions  UA, |  <  1  may  be  tested  as  a  validity 
check  for  the  inputs.)  Define  ah,  =  kh,  =  0  when  Mh  <  t  <  M.  The  results  are  two 

vectors  ah  =  (a/,o=L  ah\ . °hm)  and  kh  =  (kh],  .  .  .  ,  khM)  and  a  gain  Eh  for  each 

h. 

Step  3.  Compute  the  families  of  Lagrange  multipliers 

\0  _  A  0  \0  \  0  i 

Kh  \*h0*  A/il>  •  •  > 

corresponding  to  R°  from  ah  and  Eh  for  each  h  by  Eq.  (7)  (with  A°  in  place  of  Xh). 

Step  4.  (Initialize  variables.)  By  a  standard  LPC  algorithm,  determine  a  family 

(/?A0,  Rh\ . RhM^  of  autocorrelations  corresponding  to  the  gain  and  reflection 

coefficients  Eh  and  kh  for  each  h\  use  these  as  initial  values  for  the  variables  Rh 
(This  results  in 

Rh,  =  K 

for  0  <  t  <  Mh\  hence  computing  Rh  from  Eh  and  kh  is  actually  unnecessary  unless 
Mh  <  M.)  Also,  set 

xh  —  x° 

and 

P,  *—  0,  t  “  0,  1  ,  . . .  ,  m. 

Step  5.  (Beginning  of  main  loop.)  Compute  coefficients 

&h  =  ^  bh0=  1 ,  bh  |  ....  ,  bh  2a/) 
from  ah  for  each  h  by  Eq.  (27). 

Step  6.  Compute  the  quantities 

Sh  =  (ShO'  $h\ . 

for  each  h  by  an  algorithm  for  computing  autocorrelations  from  inverse  filter 
coefficients  and  gains.  As  inputs,  let  the  family  b„  of  coefficients  from  Step  5  play  the 
role  of  inverse  filter  coefficients  and  the  squared  power  Ef  play  the  role  of  gain.  The 
family  Sh  comprises  the  resulting  "autocorrelations."  Define  the  Jacobian  matrices  J/, 
in  terms  of  Sh  by  Eq,  (26). 
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Step  7.  Compute  the  quantities  J^lxj’  +  D0  -  X/,)  Define  A0  by  Eq.  (32),  and  replace  0 
by  a  new  value, 

0  -  0  +  A0. 

Obtain  a  tentative  value  for  ARh  for  each  h  from  Eq.  (31).  (At  this  point  the  con¬ 
straint  equation  (Eq.  (30)]  should  be  satisfied.) 

Step  8.  Check  that  the  sums  Rh  +  A Rh  are  valid  families  of  autocorrelations.  We  can  do  this 
by  using  the  Levinson  recursion  to  compute  reflection  coefficients  and  checking  that 
these  are  less  than  1  in  magnitude;  it  is  also  necessary  to  check  that  each  Rh0  is  posi¬ 
tive.  If  Rf i  +  A Rh  fails  the  test  for  any  h ,  then  reduce  the  A Rh  for  all  h  by  a  constant 
factor.  Repeat  the  testing  and  reduction  until  the  Rh  +  A Rh  pass  the  test.  Then 
replace  the  Rh  by  new  values 

Rh-  Rh  +  &Rh 

for  each  h. 

Step  9.  Decide  whether  to  continue.  If  the  new  value  of  R/,  (computed  in  Step  8)  is  not 
equal  to  the  previous  value  for  each  h  to  within  a  specified  relative  tolerance,  the 
main  loop  will  have  to  be  repeated  from  Step  5.  Even  if  the  new  Rh  are  close  enough 
to  the  previous  values,  the  main  loop  should  be  repeated  if  Rh  had  to  be  reduced  in 
Step  8.  If  the  loop  is  to  be  repeated,  continue  with  Step  10  before  going  back  to  Step 
5.  If  results  from  Step  10  are  required  for  output,  do  Step  10  before  exiting.  Other¬ 
wise  terminate  the  algorithm  now. 

Step  10.  At  this  point  new  values  for  the  Rh  from  Step  8  are  available.  Corresponding  new 
values  for  the  quantities  ah,  Eh ,  and  kh  will  be  needed  if  there  is  to  be  another  itera¬ 
tion  of  the  main  loop.  They  may  also  be  desired  as  outputs  from  the  algorithm. 
Nothing  further  need  be  done,  however,  if  there  is  not  to  be  another  iteration,  and  if 
the  only  outputs  desired  are  the  posterior  autocorrelation  estimates.  Values  for  ah 
and  Eh  may  already  be  available  as  results  of  the  use  of  the  Levinson  recursion  in 
Step  8;  if  not,  they  should  be  computed  now.  To  obtain  XA  from  ah  and  Eh ,  use  Eq. 
(7).  Then  either  go  back  to  Step  5  or  terminate  the  algorithm,  whichever  is  called 
for. 


DISCUSSION 

Both  algorithms  have  been  implemented  in  APL  and  in  FORTRAN.  The  APL  implementations 
were  written  first  for  ease  of  programming  and  debugging.  The  APL  functions  were  then  converted  to 
FORTRAN  for  incorporation  in  a  commercial  signal-processing  package  that  we  use  for  experimental 
work  with  speech.  Both  versions  of  the  algorithms  are  available  from  the  author. 

The  second  algorithm  is  so  much  faster  than  the  first,  that  the  first  can  be  recommended  only  for 
applications  that  really  demand  the  generality  of  unequally  spaced  lags,  unequally  spaced  frequencies,  or 
completely  arbitrary  prior  spectral  shapes.  The  FORTRAN  version  of  the  first  algorithm,  running  on  a 
VAX  1  1/750*  with  floating-point  accelerator,  uses  about  18.5  s  of  cpu  (cenlral-processing-unit)  time 
per  frame  of  speech  in  a  typical  speech-processing  application  —  a  two-signal  analysis  with  prior  spectral 
estimates  given  at  128  equispaced  frequencies  and  autocorrelation  constraints  given  at  13  lags: 
0,  1 . 12.  The  FORTRAN  version  of  the  second  algorithm,  running  on  the  same  machine,  uses 
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about  2.0  s  per  frame  on  the  same  problem,  where  the  prior  spectral  estimates  are  12lh-order  all-pole 
spectra  specified  by  autocorrelation  values  at  lags  0,  1 . 12. 

The  second  algorithm  also  uses  less  space  than  the  first.  In  terms  of  the  number  p  of  signals,  the 
number  n  of  discrete  frequencies,  and  the  order  m  of  the  autocorrelation  constraints,  the  FORTRAN 
implementation  of  the  first  algorithm  uses  array  storage  of  3(m  +  l)n+  (m  +  l)2+  3 pn  + 
in  +  4(/w  +  l)  floating-point  (or  double-precision)  locations,  2(w  +  l)  integer  locations,  and  an  addi¬ 
tional  amount  independent  of  the  dimensions  of  the  inputs.  For  most  applications,  the  dominant  term 
is  3 (m  + 1)«,  which  corresponds  to  storage  for  the  matrices  C,  N,  and  N'f  of  Eqs.  (2),  (18),  and  (23). 

In  terms  of  /?,  m,  and  the  order  M  of  the  all-pole  posterior  spectra,  the  second  algorithm  uses 

10(A/+l)p  +  (M+l)2  -I-  7(A/+1)  +  5(m  +  l)  floating-point  locations,  2 (m  -F  1 )  integer  locations, 
and  a  fixed  additional  amount.  Typical  values  of  the  size  parameters  are  p  =  2,  n  =  128,  and 

m  —  M  =  12.  With  these  values,  the  first  algorithm  uses  6365  floating-point  locations,  as  compared  with 

585  for  the  second.  Of  the  6365  locations,  nearly  80%  are  occupied  by  the  three  matrices  that  account 
for  the  dominant  term  3(m  +  l)n. 

We  do  not  claim  to  have  produced  optimally  efficient  multisignal  MCESA  algorithms;  there  may 
exist  algorithms  much  superior  to  the  two  presented  here  in  their  respective  domains.  The  efficiency  of 
these  algorithms  should  not,  therefore,  be  taken  as  a  measure  of  the  intrinsic  efficiency  of  multisignal 
MCESA.  The  second  algorithm,  however,  is  adequate  to  permit  considerable  experimentation  with 
speech  noise  reduction,  and  it  interfaces  well  with  conventional  LPC  analysis/ synthesis  software.  It  has 
been  possible  to  process  entire  sentences  of  speech  corrupted  with  real  or  synthetic  additive  noise,  to 
synthesize  speech  from  the  posterior  speech  autocorrelation  estimates,  and  to  listen  to  the  results. 
Effective  noise  suppression  has  been  demonstrated  in  experiments  with  a  variety  of  methods  for  choos¬ 
ing  prior  speech  and  noise  estimates  and  estimating  autocorrelations  from  the  sum  of  the  speech  and 
the  noise  signals  [3], 
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Appendix  A 

APL  PROGRAM  FOR  FIRST  ALGORITHM 


The  argument  P  of  the  APL  function  MCESP2P  below  corresponds  to  /  and  P  in  the  conven¬ 
tional  notation.  The  argument  R  corresponds  to  /  and  Rt01,  and  the  result  Q  corresponds  to  / and  Q.  P 
is  a  matrix  with  frequencies  fk  in  row  1;  each  additional  row  contains  prior  spectral  estimates  Phk  for 
one  signal  component  h.  R  is  a  two-rowed  matrix  with  lags  tr  in  row  1  and  autocorrelations  Rxox  in  row 
2.  Q  is  a  matrix  with  the  same  shape  and  format  as  P  the  frequencies  fk  in  row  1  and  posterior  spec¬ 
tral  estimates  Qhk  in  the  remaining  rows. 


MCESP2PW  is  a  version  of  the  program  that  implements  a  generalization  of  multisignal  MCESA 
that  is  described  elsewhere  [Al],  The  argument  P  has  an  additional  first  column  containing  "relative 
weights."  The  first  element  of  the  column  is  ignored.  Each  other  element  of  the  first  column  is  a 


weighting  factor  associated  with  the  prior  spectral  esti 


mate  in  the  rest  of  the  same  row.  The  format  is 


* 

A  - 

A 

w, 

/>„  ... 

Pu 

WP 

- 

Ppn 

The  argument  R  and  the  result  Q  are  as  for  MCESP2P, ;  £)has  no  additional  first  column. 


REFERENCE 


Al.  R.W.  Johnson  and  J.E.  Shore,  "Multisignal  Minimum-Cross-Entropy  Spectrum  Analysis  with 
Weighted  Priors,"  NRL  Report  8731,  in  publication;  submitted  to  IEEE  Trans.  Acoust.,  Speech. 
Signal  Process. 
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7  Q+P  MCESP2P  R i C i F i X i L i QO ; & L 

[1]  C’*-2x20(ff[l  ;]xo2)<>  .  xF-«-p[i ;  ] 

[2]  R-*-Rl  2  ;  ] 

[3]  P+-  1  0  +  P 

C  4  ]  *-/?+.  xS$C 
C  5  3  £.«-  (  p  R )  p  0 

[6]  P+iQ-P 

[7]  A-.QO+Q 

[8]  /?-(+/«*  2)*0.5 

[9]  AZ>(  (  (+/«)-*)*/?)  BS^XpOpR 

[10]  B_ :  Q+i  P+  (  pP )  p  (  L  +  A  L )  +  .  *C 

[11]  -(a/,G>0)/£ 

[12]  AZ>AL*0.9*1-L/,30*G 

[13]  -*-B_ 

[14]  CiL+L+M 

[15]  -(*/ ,Q*Q0)/A 

[16]  G*-F,[l]  Q 
V 


7  Q+P  MCESP2PV  R  i  W  i  F •,  C  i  X  i  L  i  Q0  i  h L 

[1]  V-*-*l+P[;l] 

[2]  F-»- 1  +  P  [  1  ;  ] 

[3]  P*2x20(P[ 1 ; ] xo2 ) « . xF 

[4]  F-R[2;] 

[5]  P«-  1  1  +P 

[6]  *«-/?+.  xS$P 

[7]  pR)  pO 

[8]  P+*Q+-P 

[9]  A_i  Q0+-Q 

[10]  R«-(tf+.  xQ*2)*0. 5 

[11]  AZ>(  ((+/C) -*)*/?)  1^7*  (pC)pR 

[12]  BiQ-iP+Wo.  x(Z,+  A£)  +  .  xC 

[13]  -*•(*/, Q>0)/£ 

[14]  AZ>Al,x0.9*l-L/ ,Q0*S 

[15]  *B 

[  16  ]  £:  L+-L+  AL 

[17]  -*(v/,<2*<20)/4 

[18]  Q-F,[l]  C 


Appendix  B 

FORTRAN  PROGRAM  FOR  FIRST  ALGORITHM 

The  following  FORTRAN  subroutine  DOMC2  implements  a  two-signal  version  of  the  first  algo¬ 
rithm.  It  will  handle  autocorrelations  of  order  up  to  29  and  up  to  513  discrete  frequencies.  The  array  C 
should  be  initialized  in  the  calling  program  by 

DO  1  I  =  1,M 
DO  1  K-l.N 

1  C(I,K)  =  2  *  COS (2  *  PI  *  T (I)  *  F«0) 

before  the  first  call  on  DOMC2.  It  need  not  be  reinitialized  for  subsequent  calls  as  long  as  the  same 
sets  of  frequencies  F(K)  and  lags  T (I)  are  to  be  used.  The  input  argument  WT  ("relative  weights")  is 
concerned  with  a  generalization  of  multisignal  MCESA  that  is  described  elsewhere  [B 1  ] ;  if  WT(1)  and 
WT(2)  are  made  equal  and  positive  (say  both  equal  to  1.0),  the  algorithm  will  behave  as  here 
described. 

DOMC2  requires  a  function  DISTD  and  a  subroutine  MPD,  both  shown  below.  DISTD  defines 
the  stopping  criterion. 

MPD  defines  its  argument  INV  as  the  Moore-Penrose  generalized  inverse  of  the  transpose  of  its 
argument  C,  provided  that  the  rows  of  C  are  linearly  independent.  Lawson  and  Hanson  [B2]  give 
general-inverse  routines  that  are  numerically  superior,  though  more  complex.  MPD  uses  a  matrix 
inversion  subroutine  MINVD,  not  shown.  The  form  of  the  call  is  CALL  MINVD(A,N,D,L,M),  where 
A  is  an  N-by-N  double-precision  matrix,  and  L  and  M  are  integer  scratch  vectors  of  length  N.  (A 
determinant  may  be  returned  through  D  but  is  ignored  in  MPD.)  A  contains  the  input  matrix,  which  is 
destroyed  and  replaced  by  the  inverse. 

REFERENCES 

Bl.  V.  Johnson  and  J.  E.  Shore,  "Multisignal  Minimum-Cross-Entropy  Spectrum  Analysis  with 
.ghted  Priors,"  NRL  Report  8731,  in  publication;  submitted  to  IEEE  Trans.  Acoust.,  Speech, 
Signal  Process. 

B2.  C.L.  Lawson  and  R.J.  Hanson,  Solving  Least-Squares  Problems,  Prentice-Hall,  Englewood  Cliff's, 
N.J.,  1974. 
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SUBROUTINE  DOMC2 ( C, R, P, M , N, WT , B , Q, INV, N 1 , IERR, LUN ) 

REAL*8  R(30),B(30),DB(30),P(513,2),Q(513,2),W(513),X(513) 
REAL*8  C ( M, N ) , INV (M,N) ,N1(M,N) ,TMP(30,30) , QO ( 5 1 3 , 2 ) , WT ( 2 ) 
REAL*8  AVG, D, QSUM , QSUMSQ, RMI5 13, SCALE 
INTEGER  K1 ( 30 ) ,K2 ( 30 ) 

LOGICAL  FLAG 

C  SUBROUTINE  TO  PERFORM  MINIMUM  CROSS-ENTROPY  SPECTRAL  ANALYSIS 
C  ON  MULTIPLE  SPECTRA 
C  DOUBLE  PRECISION  VERSION 
C 

C  NAVAL  RESEARCH  LABORATORY 

C  J.  T.  BUCK 

C  ADAPTED  FROM  APL  VERSION  BY  R.  JOHNSON 

C 

C  INPUT  ARGUMENTS: 

C  C  -  ARRAY  OF  COSINE  TERMS  -  2* COS ( 2*PI *T ( I ) *F( K ) ) 

C  -  FOR  TIME  AND  FREQUENCY  POINTS 

C  R  -  AUTOCORRELATIONS 

C  P  -  PRIOR  SPECTRAL  ESTIMATE 

C  FOR  EACH  OF  2  SPECTRA 

C  M  -  NUMBER  OF  TIME  POINTS 

C  N  -  NUMBER  OF  FREQUENCY  POINTS 

C  WT-  RELATIVE  WEIGHTS  OF  PRIOR  SPECTRA 

C  LUN  IF  >  0,  LUN  FOR  STATUS  MESSAGES 

C  OUTPUT  ARGUMENTS: 

C  B  -  LAGRANGE  MULTIPLIERS  USED  IN  SOLUTION 

Q  -  MINIMUM  CROSS-ENTROPY  SPECTRAL  ESTIMATES  FOR  EACH  SPECTRUM 
IERR  -  ERROR  FLAG  -  0  IF  OK,  NONZERO  OTHERWISE 

SCRATCH  ARGUMENTS: 

N 1  -  WORK  MATRIX 

INV  -  GETS  GENERALIZED  INVERSE  OF  WORK  MATRIX 
W,X  -  WORK  VECTORS 
C  QO  -  SAVES  PREVIOUS  VALUE  OF  Q 

C  DB  -  INCREMENT  FOR  CALCULATING  NEXT  B 

C  IMP , K1 , K2  -  TEMPORARIES  REQUIRED  FOR  MPD  ROUTINE 

C 

C  SUBROUTINES  REQUIRED 

C  MPD  -  CALCULATES  GENERALIZED  INVERSE 

C  DISTD  -  DISTANCE  MEASURE 

C  MINVD  -  INVERTS  A  MATRIX  (CALLED  BY  MPD) 

C 

IERR=0 

C  GET  AVG  VALUE  OF  PRIORS  (FOR  CONVERGENCE  CHECK) 

AVG=0 . 0 
DO  1=1, N 
DO  J=1 , 2 

AVG=AVG+P ( I , J ) 

ENDDO 

ENDDO 

AVG=AVG/(N*2) 

C  INITIALIZE  B , C , AND  X 
DO  J=1 ,M 
B(  J  )=0 . 0 
ENDDO 
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CALL  MPD( C , INV, M, N, TMP, K1 , K2 ) 

DO  K=1,N 
DO  L=1 , 2 

Q( K, L)=P (K, L ) 

P(K,L)=1. 0/P (K, L) 

ENDDO 
X(K)=0.0 
DO  J=1 , M 

X(K)=X(K)+INV(J,K)*R(J) 

ENDDO 
ENDDO 
NPASS=0 
NRSOJ=0 

MAIN  LOOP  STARTS  HERE 
NPASS=NPASS+1 
IF ( NPASS . GT . 100 )THEN 

C  Convergence  failure  -  type  message,  return  error 
D=DISTD (Q,Q0,N*2 ) /AVG 
TYPE  3 1 , D 

31  FORMAT ( '  Convergence  failed  after  100  passes  -  D=  '  ,1PE13.6) 

IERR= 1 
RETURN 
ENDIF 

C 

DO  K=  1 ,  N 

C  COPY  OLD  VALUE  OF  Q,  CALC  AVG,  MEAN  SQUARE  VALUES 
QSUM=0 . 0 
QSUMSQ=0 . 0 
DO  L=1 , 2 

Q0(K,L)=Q(K,L) 

QSUM=QSUM+Q(K,L) 

QSUMSQ=QSUMSQ+Q (K,L) *  *2/WT { L ) 

ENDDO 

QSUMSQ=SQRT ( QS  UMSQ ) 

W( K )= (QSUM-X ( K ) J/QSUMSQ 
DO  J=1 ,M 

N1(J,K)=C(J,K)*QSUMSQ 

ENDDO 

ENDDO 

C  GET  INVERSE  OF  WORK  MATRIX 

CALL  MPD ( N 1 , INV , M , N , TMP , K 1 , K2 ) 

C  CALCULATE  INCREMENT  FOR  LAGRANGE  MULTIPLIERS 
DO  J=1,M 
DB( J)=0.0 
DO  K=  1 ,  N 

DB(J)=DB(J)+INV(J,K)*W(K) 

ENDDO 

ENDDO 

C  Calculate  new  Q.  Check  for  overshoot.  Loop  until  no  overshoot  occurs. 
FLA G=. TRUE. 

DO  WHILE (FLAG) 

FLAG*. FALSE. 

C  Calculate  new  Q's  from  P's  and  B’s.  If  any  are  negative,  set  FLAG 
DO  K= 1 , N 
DO  L=  1 , 2 
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Q(K,L)=P(K,L) 

DO  J=  1 , M 

Q(K,L)=Q(K,L)+(B(J)+DB(J) )*C( J,K)/WT(L) 

ENDDO 

Q(K,L)=1.0/Q(K,L) 

IF ( Q (K , L ) . LT . 0 . 0 ) FLAG= . TRUE . 

ENDDO 

ENDDO 

C  If  overshoot  occured  (resulting  in  negative  sprectral  power) 
IF  (FLAG)  THEN 

C  Then  scale  the  increment  and  repeat 
RMIN=1 . E37 
DO  K=  1 ,  N 
DO  L= 1 <  2 

IF(Q0(K,L)/Q(K,L) . LT. RMIN ) RMIN=Q0 (K, L) /Q ( K, L ) 
ENDDO 
ENDDO 

SCALE=0 . 9/ ( 1 .0-RMIN) 

NRSOJ=NRSOJ+1 
DO  J=1 ,  M 

DB( J )=DB ( J )* SCALE 
ENDDO 
ENDIF 
ENDDO 

C  Calculate  new  b  value ,test  for  termination 
DO  J= 1 , M 

B  ( J ) =B ( J ) +DB ( J ) 

ENDDO 

IF(DISTD(Q,Q0,N*2).LT.1.E~10*AVG)THEN 
IF ( LUN . NE . 0 ) WRITE ( LUN ,100) NPASS , NRSOJ 
RETURN 
ENDIF 

100  FORMAT( 12X, 'Passes : ‘ , 13, '  Overshoots 13 ) 

GOTO  30 
END 


FUNCTION  DISTD( A, B, N) 

REAL*8  A( N ) , B( N ) , DISTD 
C 

C  J.  T.  BUCK 

NAVAL  RESEARCH  LABORATORY 

C  THIS  ROUTINE  CALCULATES  A  DISTANCE  MEASURE  BETWEEN  THE  VECTORS 
C  A  AND  B,  WHICH  IS  SIMPLY  THE  MAXIMUM  ABSOLUTE  VALUE  DIFFERENCE 
C  OF  THE  COMPONENTS. 

DISTD=-1 .E37 
DO  10  1=1, N 

IF ( ABS ( A( I )-B( I ) ) .GT. DISTD )DISTD=ABS(A( I ) -B ( I ) ) 

10  CONTINUE 

RETURN 
END 
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SUBROUTINE  MPD ( C, INV,M, N, TMP, K1 , K2 ) 

C 

C  MPD:  EMULATES  APL  DOMINO  OPERATOR  (GENERALIZED  INVERSE) 

C  DOUBLE  PRECISION  VERSION 

C 

C  J.  T.  BUCK 

C  NAVAL  RESEARCH  LABORATORY 

C 

C  INPUT  ARGUMENTS: 

C 

C  C  -  M  BY  N  MATRIX  TO  INVERT 

C  M,  N  -  MATRIX  DIMENSIONS 

C 

C  OUTPUT  ARGUMENTS: 

C 

C  INV  -  M  BY  N  RESULT 

C 

C  SCRATCH  ARGUMENTS: 

C 

C  TMP  -  M  BY  M  REAL  MATRIX 

C  K1,K2  -  INTEGER  INDEX  VECTORS,  LENGTH  M 

C 

C  SUBROUTINES  REQUIRED: 

C 

C  MINVD  -  INVERTS  A  MATRIX 

C 

REAL*8  C ( M, N ) , INV( M, N) , TMP( M,M) , D 
INTEGER  K1(M),K2(M) 

C 

C  THIS  ROUTINE  CALCULATES  THE  GENERALIZED  INVERSE  OF  THE  TRANSPOSE 
C  OF  C  AND  RETURNS  THE  RESULT  IN  INV.  THE  ARGUMENTS  TMP, K1, AND  K2 
C  ARE  USED  BY  MINVD,  WHICH  INVERTS  A  MATRIX. 

C 

C  INV=INVERSE ( C  TIMES  C-TRANSPOSE )  TIMES  C 
C 

C  CALCULATE  TMP=C  TIMES  C-TRANSPOSE 
DO  10  1=1, M 
DO  10  J=1 ,M 
TMP(I,J)  =  0.0 
DO  10  K=1,N 

TMP ( I , J ) =TMP ( I , J ) +C (I,K)*C(J,K) 

1 0  CONTINUE 

C  GET  THE  INVERSE  OF  IMP 

CALL  MINVD (TMP, M, D, K1 , K2 ) 

C  CALCULATE  INV=  TMP  TIMES  C 
DO  20  1=1, M 
DO  20  J=  1 ,  N 
INV( I , J )=0 . 0 
DO  20  K=  1 ,  M 

INV( I , J)=INV( I , J)+TMP( I , K)*C( K, J) 

20  CONTINUE 

RETURN 
END 
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APL  PROGRAM  FOR  SECOND  ALGORITHM 


The  arguments  RO  and  RTOT  and  the  result  R  of  the  APL  function  MCEAUTW  below 
correspond  to  the  R®,  A'01,  and  the  Rh  in  conventional  notation.  RO  is  a  matrix  with  one  row  for  each 
signal  component  h  and  one  column  for  each  autocorrelation  lag  from  0  through  Mh  (which  is  indepen¬ 
dent  of  h\  the  APL  version  requires  all  the  priors  to  be  of  the  same  order).  In  addition,  column  1  con¬ 
tains  "relative  weights,"  which  are  concerned  with  a  generalization  of  multisignal  MCESA  that  is 
described  elsewhere  [Cl].  Thus  the  prior  autocorrelation  information  occupies  columns  2  through 
Mh  +  2  of  RO.  The  user  who  is  not  concerned  with  relative  weights  may  simply  set  all  elements  of  the 
first  column  of  RO  equal  to  1.  RTOT  is  a  vector  of  length  m  +  1.  R  is  returned  as  a  matrix  with  the 
same  number  of  rows  as  RO  and  one  column  for  each  lag  from  0  through  M  (but  no  column  of 
weights).  MCEAUTW  uses  several  other  functions  for  discrete  convolution  and  various  LPC  parameter 
conversions. 

If  the  arguments  A"  and  K  of  CVR  are  vectors  of  length  M  and  length  N,  respectively,  the  result  is 
their  discrete  convolution,  a  vector  of  length  M+N—  1.  If  X  and  Y  are  matrices  with  the  same  number 
of  rows,  the  result  is  a  matrix,  also  with  the  same  number  of  rows;  each  row  of  the  result  contains  the 
discrete  convolution  of  the  corresponding  row  of  X  with  the  corresponding  row  of  Y.  CVR  similarly 
extends  row  by  row  to  arrays  with  any  number  of  dimensions.  Thus  CVR  applied  to  a  2x3x4  array  and 
a  2x3x5  array  yields  a  2x3x8  result. 

The  functions  for  LPC  parameter  conversions  likewise  extend  systematically  to  higher¬ 
dimensional  arguments,  although  we  describe  only  their  application  to  arguments  of  the  lowest  dimen¬ 
sion.  Thus  AUTAPAR  is  said  to  apply  to  a  vector  K  and  a  scalar  RO  to  yield  a  vector  result  R.  How¬ 
ever,  it  is  then  to  be  understood  that  AUTAPAR  may  also  apply  to  a  matrix  K  and  a  vector  RO  with 
one  element  for  each  row  of  K  to  yield  a  matrix  result  R  with  one  row  for  each  row  of  K. 

The  argument  to  PREAPAR  is  a  vector  K  containing  parcor  coefficients,  or  negative  reflection 
coefficients,  The  result  PREAPAR  K  is  a  vector  A  containing  inverse  filter 

coefficients.  The  coefficient  a0  is  absent,  and  the  sign  convention  is  the  opposite  of  that  used  in  the 
FORTRAN  programs  and  the  main  body  of  this  report.  Thus  the  contents  of  A  are  -a]t  .  .  .  ,  ~aM. 

The  function  PARAPRE  is  an  inverse  to  PREAPAR.  The  argument  is  a  vector  A  of  inverse  filter 
coefficients  (with  reversed  sign),  and  the  result  PREAPAR  A  is  a  vector  K  of  negative  reflection 
coefficients. 

The  function  LAGAPRE  takes  as  arguments  a  vector  A  of  inverse  filter  coefficients  (with  reversed 
sign)  and  a  scalar  gain  £  The  result  £  LAGAPRE  A  is  computed  by  Eq.  (7),  but  without  the  factor  of 
'h  in  \0  (and  we  have  suppressed  the  subscript  h).  The  result  is  thus  a  vector  B  containing 
2*o.  *i . *w- 

The  arguments  to  AUTAPAR  are  a  vector  K  containing  negative  reflection  coefficients  and  a 
scalar  RO  containing  a  total  power  R0.  The  result  RO  AUTAPAR  K  is  a  vector  R  containing  autocorre¬ 
lation  values  R0,  R\,  .  .  .  ,  R w. 

PARAAUT  an  inverse  function  to  AUTAPAR ,  takes  R  as  its  argument  and  returns  K  as  its  result 
PARAAUT  R. 


RODNEY  W.  JOHNSON 


GAISPAR  yields  a  scalar  result,  the  ratio  of  LPC  gain  to  total  power,  or  E/R 0  in  the  notation  of 
Eq.  (25)  (with  subscripts  h  suppressed).  This  result  is  computed  as  a  function  of  negative  reflection 
coefficients  K. 

REFERENCE 
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Weighted  Priors,"  NRL  Report  8731,  in  publication;  submitted  to  IEEE  Trans.  Acoust.,  Speech. 
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7  R*-R 0  MCEAUTW  RTOT; W \ K\D \MU \M\F \ II \I2 \ A\E\ LO \L\B\ S\T \LR \J \ EB\ SC ALED 

[1]  W+-t  RO  [  ;  1  ] 

[2]  B0«-  0  1  +i?0 

[3]  K+PAR  EAUT  RO 

[4]  -••(  lv.5|  (  ,K)  .PARbAUT  RTOT)/0 

[5]  £>*-(  pflO  )  [  1  ] 

[6]  MU-'l+pRTOT 

[7]  1  + (pfl0)[2] 

[8]  F+-(D  ,M+  l)p0.5.«pl 

[9]  Il+0,iMU 

[10]  I2-*-(  1 1+M+l)  °  - 1 1 

[11]  Il«-(  I 1+M+l ) ° . +11 

[12]  A-PRE A PAR  K 

[13]  £>ff0[  il]*GAI t>PAR  K 

[14]  L0+-F*(.D  ,M+1  )+£  LAG  hPRE  A 

[15]  A--{DtM)*A 

[16]  L+-L0 

[17]  B*-(.M+1  )p0 

[18]  R+-R0 

[19]  -►(  (p/?0)[2]2W+l  )/A 

[20]  R+RO  [  ;  1  ]  AUTEPARU)  ,M)iK 

[21]  A  iK+-PARbPRE  0  1  +-(~1,/1)  CVR  "l  ,A 

[  22]  S«-(  (E*2)*GAILPAR  K)  AUTLPAR  K 

[23]  S«-(4>  0  1  +(0  ,-M)*S)  ,S 

[24]  T+-  (  0  ,  Af )  +  (  0 ,-M) IS  CTfl4>U0+Vo.xBW, 

[25]  A/?^-(4>(0,-«)*D  +  (0,#)*r 

[26]  J*--V+.*S 

[27]  EB+-{RT0T-(,MU+1  )  WB+A/f)  S</[Il  W[I2] 

[28]  B*-B+  (M  +  l  )  +  AB 

[29]  T+(.0,MU)*(0,-M)+S  CVR  V°.*<|>AB 

[30]  A/?^-A/?-(4>(0,-W)+r)  +  (0,M)+r 

[31]  5C/1BE0-O 

[32]  B0«-fl 

[33]  -*•(([  /  \RTOT-  +  JRO  +  bR)<,RTOTll]*lE-5)/B 

[34]  » CONSTRAINTS  VIOLATED' 

[35]  -*-0 

[36]  B  :  B-*-i?0  +  AB 

[37]  K+-PAR tiAUT  R 

[38]  -*■(  (  Oa.  <S[  ;1  ])aia.>  |  ,K)/C_ 

[39]  AB-*-Aff  x  0  .  7 5 

[40]  SCALED+1 

[41]  -*B 

[42]  C:*(SCALED)/D 

[43]  A/ ,R=R0 )/0 

[44]  D:A+PRELPAR  K 

[45]  E+Rl  HlxGAlLPAR  K 

[46]  Z>Fx£  LAGLPRE  A 

[47]  -+A 
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7  Z+Y  CVR  XiNiNiR 

Cl] 

M+~ lipX 

C2] 

N+~lipY 

• 

[3] 

fl-ppX 

[4] 
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Appendix  D 

FORTRAN  PROGRAM  FOR  SECOND  ALGORITHM 


The  following  FORTRAN  subroutine  MCLPC  implements  a  one-  or  two-signal  version  of  the 
second  algorithm  and  will  handle  autocorrelations  of  order  29  or  less.  The  input  argument  W  (relative 
weights)  is  concerned  with  a  generalization  of  multisignal  MCESA  that  is  described  elsewhere  [Dl],  if 
W(l)  and  W(2)  are  made  equal  and  positive  (say  both  equal  to  1.0),  the  algorithm  will  behave  as  here 
described. 

MCLPC  requires  a  matrix-inversion  subroutine  MINVD  and  subroutines  DLEVIN,  DA2L, 
DRC2AU,  and  DA2RC  for  conversions  among  various  sets  of  LPC  parameters. 

The  form  of  a  call  for  MINVD  is  CALL  MINVD  (A,N,D,L,M),  where  A  is  an  N-by-N  double¬ 
precision  matrix  and  L  and  M  are  integer  scratch  vectors  of  length  N.  (A  determinant  may  be  returned 
through  D  but  is  ignored  in  MCLPC.)  A  contains  the  input  matrix,  which  is  destroyed  and  replaced  by 
the  inverse. 

DLEVIN  implements  the  standard  Levinson  recursion.  A  call  for  DLEVIN  has  the  form  CALL 
DLEVIN  (MP,R,A,ALPHA,RC),  where  R,  A,  and  RC  are  double-precision  vectors,  ALPHA  is  a 
double-precision  scalar,  and  MP  is  an  integer.  The  inputs  are  MP,  which  is  the  order  plus  one,  and  R, 
which  contains  MP  autocorrelation  values  from  lag  0  to  lag  M  =  MP— 1.  The  outputs  are  MP  inverse 
filter  coefficients  in  A,  M  reflection  coefficients  in  RC,  and  the  gain  (excitation  power)  in  ALPHA. 
A ( 1 )  contains  a0=  1.  The  subroutine  body  is  adapted  from  Market  and  Gray’s  subroutine  AUTO  [D2, 
p.  51]  with  the  following  changes:  convert  to  double  precision;  redimension  local  arrays  to  60  (the  sub¬ 
routines  must  handle  twice  the  order  of  the  input  autocorrelations);  change  MP  =  M-f-l  to  M  =  MP— 1; 
remove  the  initial  nested  DO  loops  that  compute  R  from  signal  samples. 

DA2L,  shown  below,  computes  Lagrange  multipliers.  The  inputs  are  MP1,  the  order  plus  1;  an 
array  A  of  inverse  filter  coefficients;  and  a  gain  (excitation  power)  EIN.  The  output  XL  contains  MP1 
Lagrange  multipliers  computed  by  Eq.  (7)  except  for  the  factor  of  ]h  in  kh0  =  XL(1). 

DRC2AU,  shown  below,  computes  autocorrelations  from  reflection  coefficients  and  the  square 
root  of  the  gain.  DRC2AU  uses  a  subroutine  DRC2E  for  computing  the  total  power  E  from  reflection 
coefficients  and  the  square  root  of  the  gain. 

DA2RC  computes  reflection  coefficients  from  inverse  filter  coefficients.  A  call  has  the  form 
CALL  DA2RC(A,RC,M,IERR),  where  A  and  RC  are  double-precision  vectors,  and  M  and  IERR  are 
integers.  The  inputs  are  the  order  M  and  M  +  l  inverse  filter  coefficients  in  A.  The  outputs  are  M 
reflection  coefficients  in  RC  and  an  error  indication  in  IERR;  a  nonzero  value  for  IERR  indicates  that  a 
computed  magnitude  for  a  reflection  coefficient  was  not  less  than  1.  The  subroutine  body  is  adapted 
from  Markel  and  Gray’s  subroutine  STEPDN  [D2,  p.  96]  with  the  following  changes:  convert  to  dou¬ 
ble  precision  and  redimension  local  arrays  as  for  DLEVIN,  and  return  the  error  indication  through 
IERR  instead  of  printing  a  message. 
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SUBROUTINE  MCLPC( RIN, R,W/MDP1 , NS, MPR, A, RC , E , L, JAC , IERR ) 

C  ■  •  < 

C...  Multiple  signal  minimum  cross-entropy  spectral  analysis 

C...  with  LPC  priors 

C*  •  t 

C . . .  J.  Buck 

C...  Based  on  APL  program  by  R.  Johnson 

C  •  •  • 

C...  Double  precision  version 

C...  Includes  degree  of  belief  parameters 

C...  Allows  different  model  orders  for  priors,  input 

C «  .  • 

C...  RIN  -  Input  autocorrelations  for  total  signal 

C...  R  -  Set  of  prior  autocorrelations  -  destroyed  and  replaced 

C...  by  posterior  autocorrelations 

C...  W  -  Relative  weights  for  each  prior  estimate  (degree  of  belief) 

C...  MDP1  -  Total  number  of  input  autocorrelations  (MD+1) 

C...  NS  -  Number  of  signals 

C...  MPR  -  LPC  model  order  for  each  prior  (vector,  length  NS) 

C...  A  -  Autoregressive  coefficients  for  posterior  signals 

C...  RC  -  Reflection  coefficients  for  posterior  signals 

C...  E  -  LPC  error  powers  for  posterior  signals 

C...  L  -  Lagrange  multipliers  for  posterior  signals 

C...  JAC  -  Scratch  array  -  MDP1  by  MDP1 

C...  IERR  -  error  flag  -  nonzero  if  input  autos  bad 

C 

INTEGER  K1(30) ,K2(30) , MPR ( 2 ) 

REAL*8  RO(30,2),R(30,2),DR(30,2),A(30,2),L(30,2),LO(30,2) 

REAL*8  RC  (  2  9 , 2  )  ,  W  (  2  )' 

REAL*8  B ( 30 ) , DB ( 30 ) , RD( 30 ) , DL( 30 ) , RTOT ( 30 ) , RIN (30) 

REAL*8  JAC ( MDP 1 , MDP 1 ) ,E<2) , SS ( 60 ) , ASQ< 60 ) ,RC2(60) 

REAL*8  S (60 , 2 ) , SCLU, SCLD, D, DIFF 
LOGICAL  SCALED 
C 

C  Initialization 
C 

IERR=0 
SCLU=RIN ( 1 ) 

SCLD=1 ./SCLU 
C 

C  Get  maximum  model  order  M 
C 

MD=MDP 1 - 1 
M=MD 

DO  N=  1 ,  NS 

M=MAX ( MD , MPR ( N ) ) 

ENDDO 

MP1=M+1 

C 

DO  1=1, MP1 

RTOT ( I )=RIN(I )*SCLD 
DO  N=  1  ,  NS 

R(I,N)=R(I,N ) *SCLD 
ENDDO 
ENDDO 
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IPASS=0 

C 

C  Get  A  coefficients,  reflection  coefficients,  1 agrange  multipliers 
C  for  each  prior. 

C 

DO  N=1 , NS 

CALL  DLEVIN ( MPR ( N ) + 1 , R ( 1,N),A(1,N),E(N),RC(1,N)) 

DO  1=1 ,MPR(N ) 

IF ( ABS ( RC ( I , N ) ) .GE. 1 . 0 )THEN 

TYPE  ILLEGAL  PRIOR  AUTOCORRELATIONS ' 

IERR= 1 
RETURN 
ENDIF 
ENDDO 

C  Extend  priors  to  order  M 

IF(MPR(N) .LT.M)  THEN 
DO  I=MPR(N )  + 1  ,  M 
RC  ( I , N )=0 . 

A( 1+1 ,N)=0. 

ENDDO 

CALL  DRC2AU(SQRT(E(N ) ) , RC( 1 , N ) ,M, R( 1 , N ) ) 

ENDIF 

CALL  DA2L(A(  1 ,N) ,L( 1 ,N) ,MP1,E£N) ) 

L( 1 ,N)=0. 5*L ( 1 , N ) 

ENDDO 

C  Zero  work  vectors 
DO  1=1 ,MP1 
B(I)=0.0 
DO  N= 1 , NS 

RO ( I , N)=R( I , N ) 

LO ( I , N )=L ( I , N ) 

ENDDO 

ENDDO 

C 

C  Main  loop 
C 

20  CONTINUE 

IPASS=IPASS+ 1 

IF( IPASS.GT. 100 )  THEN 

TYPE  *, 'Convergence  failed  after  100  passes’ 

RETURN 

ENDIF 

DO  N=1,NS 

DO  1=1, MP1 

ASQ(I )=0. 

IM=2  *MP 1 -I 
ASQ( IM)=0. 

DO  K=  1 , 1 

ASQ< I )=ASQ( I )+A(K,N)*A( I-K+1 ,N) 

IF(  I.LT.MP1 )ASQ(IM)=ASQ( IM)+A ( MP1-K+ 1 , N )*A( MP 1-I+K, N ) 
ENDDO 
ENDDO 

CALL  DA2RC( ASQ, RC2 , 2*M, IERR ) 

CALL  DRC2AU( E (N),RC2,2*M,S( 1 ,  N  )  ) 
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DO  K= 1 , MP 1 

DL(K)=LO<K,N)+B(K)/W(N)-L<K,N) 

ENDDO 

DO  J=1 , MP 1 

DR(J,N)=0. 

DO  K= 1 , MP 1 

DR(J,N)=DR(J,N)-(S(J+K-1 , N )+S ( ABS ( J-K )+ 1 , N ) )*DL(K) 
ENDDO 
ENDDO 
ENDDO 

DO  1=1 , MP1+M 
SS(I)=0. 

DO  N=  1 ,  NS 

SS(I)=SS(I)-S(I,N)/W(N) 

ENDDO 

ENDDO 

DO  J= 1 , MDP 1 

DO  K= 1 , MDP 1 

JAC(J,K)=SS(J+K-1 )  +  SS  ( ABS ( J-K )+ 1 ) 

ENDDO 

ENDDO 

C 

C  Invert  Jacobian  in  place 
C 

CALL  MINVD ( JAC , MDP1 ,  D,  K1 ,  K2  ) 

SCALED  =  .FALSE. 

DO  1=1 ,MDP1 

RD( I )=RTOT( I ) 

DO  N= 1 , NS 

RD(I)=RD(I)-R(I,N)-DR(I,N) 

ENDDO 

ENDDO 

DO  1= 1 , MDP1 
DB ( I )=0 . 

DO  J= 1 , MDP 1 

DB(I)=DB(I)+RD(J)*JAC(J,I) 

ENDDO 

B( I )=B( I )+DB( I ) 

ENDDO 

C 

C  Compute  delta-R 
C 

DO  N= 1 , NS 

DO  J=1 , MP 1 

DO  K= 1 , MP 1 

DR( J ,N)=DR(J,N)-(S(J  +K- 1 , N )+S ( ABS (J-K )+ 1 , N ) ) *DB ( K ) /W( N ) 
ENDDO 
ENDDO 
ENDDO 
C 

Check  constraints 

D  =  0.0 
DO  1=1 ,  MDP1 

DIFF  =  RTOT(I) 
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DO  N=  1  ,  NS 

DIFF  =  DIFF  -  R ( I , N )  -  DR(I,N) 

ENDDO 

D  =  MAX ( D , ABS (DIFF) ) 

ENDDO 

IF  (D.GE. 1.0E-5)  THEN 
TYPE  70,  IPASS,  D 

70  FORMAT!'  Constraints  violated  on  pass' ,14, 

+  Rel.  discrepancy' , 1PE13. 6, 1 . 

RETURN 
END  IF 

Save  old  value  of  R 

DO  1=1, MP1 

DO  N=  1 ,  NS 

R0  ( I ,  N  )  =  R  ( I ,  N ) 

ENDDO 
ENDDO 

Calculate  new  R's,  test  for  legality 

20  DO  1=1, MP1 

DO  N=  1  ,  NS 

R(I,N)=R0(I, N )+DR ( I , N ) 

ENDDO 
ENDDO 
DO  N=1 , NS 

IF(R( 1 ,N) . LE . 0 . 0 ) GOTO  150 

CALL  DLEVIN ( MP 1 , R ( 1 , N ) , A( 1 , N ) , E(N ) , RC( 1 , N ) ) 

DO  1=1, M 

IF( ABS ( RC( I , N) ) .GE . 1 . 0 )GOTO  150 
ENDDO 
ENDDO 
GOTO  160 

New  autocorrelations  not  feasible;  reduce  size  of  jump 

iO  DO  1=1,  MP1 

DO  N=  1 ,  NS 

DR(I,N)=0.75*DR(I,N) 

ENDDO 
ENDDO 

SCALED  =  .TRUE. 

GOTO  120 

Get  new  Lagrange  multipliers 
.0  DO  N=  1 ,  NS 

CALL  DA2L{ A( 1 , N) , L( 1 , N ) ,MP1 , E(N ) ) 

L ( 1 ,N)=0.5*L( 1 ,N) 

ENDDO 

Convergence  check 
C 
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IF  (SCALED)  GOTO  20 
DO  1=1, MP1 

DO  N=  1 ,  NS 

IF(R0(I,N) .EQ.O. )GOTO  20 

IF ( ABS ( DR ( I , N ) /R0 ( I , N ) ) . GT . 1 . E-5 )GOTO  20 
ENDDO 
ENDDO 
DO  N=  1 ,  NS 

E(N)=E(N)*SCLU 
DO  1=1, MP1 

R(I,N)=R(I,N)*SCLU 

ENDDO 

ENDDO 

IERR=0 

RETURN 

END 


SUBROUTINE  DA2L( A, XL, MP 1 , EIN ) 
REAL*8  A ( MP 1 ) , XL ( MP 1 ) , EIN , EO 
EO=1 ./EIN 
DO  10  J= 1 , MP 1 
XL ( J )=0 . 0 
DO  5  1=1 , MP1-J+1 

XL( J )=XL( J )+A( I )*A( I+J-1 ) 
5  CONTINUE 

XL ( J )=XL( J ) *EO 
10  CONTINUE 

RETURN 
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SUBROUTINE  DRC2AU( RMSS ,  RC ,  N ,  R ) 

C.  .  . 

C...  CONVERT  GAIN  AND  REFLECTION  COEFFICIENTS  TO  AUTOCORRELATIONS 

C...  DOUBLE  PRECISION  VERSION 

C.  .  . 

C...  INPUTS: 

C  •  •  • 

C...  RMSS  -  SQRT  OF  LPC  GAIN 

C...  RC  -  REFLECTION  COEFFICIENT  ARRAY 

c***  N  -  NUMBER  OF  REFLECTION  COEFFICIENTS 

C*  •  * 

C.«..  OUTPUTS: 

C  •  •  • 

C...  R  -  ARRAY  OF  N+1  AUTOCORRELATIONS 

c  •  •  • 

REAL*8  RC( 1 ) ,R( 1 ) ,A(60) ,A1 (60) , RMSS , E 
CALL  DRC2E(RC,N, RMSS, R( 1 ) ) 

E=R( 1 ) 

DO  50  M= 1 , N 

R(M+1 )=-E*RC ( M) 

E=E  *  (  1  -RC ( M ) *  *  2 ) 

A 1 (M)  =  -RC(  M) 

IF(M.EQ. 1 )GOTO  20 
DO  10  J=1 ,  M- 1 

R(M+1 )=R(M+1 )+A( J)*R(M-J+1 ) 

10  A1(J)=A(J)+RC(M)*A(M-J) 

20  DO  40  J=1,M 

40  A( J )=A 1 ( J ) 

50  CONTINUE 

RETURN 
END 


SUBROUTINE  DRC2E ( RC , M, RMSS , E) 
REAL*8  RC (30), RMSS , E , CF 
CF=1 . 

IF ( M. EQ. 0 )GOTO  20 
DO  10  1=1, M 

10  CF=CF*(1.-RC(I)**2) 

20  E=RMSS**2/CF 

RETURN 
END 
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